{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1b2eb813",
   "metadata": {},
   "source": [
    "### Note:\n",
    "In this directory, there are two examples using PDFs: **Example 3 - Plot PDF for a station**, and **Example 4 - Calculate PDFs from PSDs**.  These two examples are provided in order to highlight different ways that you can use the PDF and PSD values that ISPAQ generates.  \n",
    "\n",
    "To be specific, the difference between the two examples are:  \n",
    "\n",
    "#### Example 3 - Plot PDF for a station:  \n",
    "Example 3 uses PDFs that already exist in the ISPAQ example database. This means that they have been calculated using an ispaq.py command with the `--output db --db_name ispaq_example.db` options.  \n",
    "\n",
    "This is a great way to do it, especially if you plan to run the PSDs and PDFs at the same time, say on some sort of regular schedule.  In that case, you might as well calculate both in the same command and store them both in the ISPAQ database for later retrieval.  \n",
    "\n",
    "Additionally, we have tried to make it simple to calculate PDFs in ISPAQ for cases where you already have PSDs for the time span you are interested in.  For example, PDFs calculation does not require seismic data since it instead reads in existing PSDs.  That means that if you, the user, have been calculating daily PSDs for the past year, you don’t need to load a year’s worth of data to calculate a year-long PDF  - you can just use the existing PSDs! By calculating that year-long PDF using ISPAQ, it will be saved to either the database or the csv file and you will be able to retrieve it later.  \n",
    "\n",
    "#### Example 4 - Calculate PDFs from PSDs:  \n",
    "Example 4 will calculate PDFs _on the fly_, meaning that they do not need to exist in the ISPAQ metric database, nor will they be saved to the ISPAQ metric database.   \n",
    "\n",
    "Why would you want to do this if you can simply use an ispaq.py command to calculate and save the PDFs in the database?  Here are a couple possible reasons: \n",
    "1) You may want to calculate PDFs on an arbitrary timeframe but don't feel the need to save the PDF values, say if you are just poking around at or investigating changes in the data and don't want to clutter the database.   \n",
    "\n",
    "2) To prevent the ISPAQ database from growing too complicated, the pdf table in the ISPAQ database is very simple and PDFs values are stored with the start and end times used to calculate that particular PDF.  If you calculate daily PDFs for a week and then additionally calculate a week-long PDF, the database will store 8 PDFs - one for each day in the week, and one that spans the entire week. This means that, even if you have used ISPAQ to calculate your arbitrary time frame, you must know the specific start and end times of the PDF that you are looking to retrieve. If you look for a time range using less-than and greater-than (<>) instead of equals-to (==) then you risk retrieving multiple PDFs, including ones that you did not intend. By using this on-the-fly method, you bypass this risk since PSDs are stored by the individual PSD (usually an hour span, can vary depending on the sample rate of the data), and only those PSDs that are needed to calculate the PDF are retrieved. \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "*Both methods are valid and can be useful in different situations.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "constitutional-detective",
   "metadata": {},
   "source": [
    "## Example 3 - Plot PDF for a station\n",
    "The intent of this series of Jupyter Notebooks is to demonstrate how metrics can be retrieved from the ISPAQ example sqlite database and provide some ideas on how to  use or plot those metrics.  \n",
    "\n",
    "This example creates a PDF plot for a station using existing ISPAQ PDF values. It requires that we have the PDF values already calculated for the target for the requested\n",
    "days, and those values should live in the ISPAQ example database. To generate PDFs, corrected PSD values must already exist. If they do not yet exist, then you can run them via (this will take several minutes):\n",
    "\n",
    "    ./run_ispaq.py -M psd_corrected -S ANMO --starttime 2020-10-01 --endtime 2020-10-16 --output db --db_name ispaq_example.db\n",
    "\n",
    "To calculate PDF values:\n",
    "\n",
    "    ./run_ispaq.py -M pdf -S ANMO --starttime 2020-10-01 --endtime 2020-10-16 --output db --db_name ispaq_example.db --pdf_interval aggregated\n",
    "\n",
    "Note: The above command will also create a PDF plot if the `pdf_type` parameter is set to 'plot' in the preference file or on the command line. The plot created in this Jupyter notebook has a different color scheme from the default plots and does not include the noise model or the max/mode/min curves.\n",
    "\n",
    "Or to calculate both PSDs and PDFs at the same time:\n",
    "\n",
    "    ./run_ispaq.py -M psd_corrected,pdf -S ANMO --starttime 2020-10-01 --endtime 2020-10-16 --output db --db_name ispaq_example.db --pdf_interval aggregated\n",
    "\n",
    "This example will assume that the above command has already been run and the PDFs already exist in the database.\n",
    "\n",
    "To begin, we need to import the necessary modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "spanish-charm",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sqlite3\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.dates import DateFormatter\n",
    "import matplotlib.dates as mdates\n",
    "import numpy as np\n",
    "import datetime"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "simple-tobago",
   "metadata": {},
   "source": [
    "Because PDFs are calculated for set frequency bins, which depend on the sample rate of the data, we create a simple function that will help us with placing our tick marks in the right location in the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "certain-stationery",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_nearest(array, value):\n",
    "    array = np.asarray(array)\n",
    "    idx = (np.abs(array - value)).argmin()\n",
    "    return idx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "coordinate-immigration",
   "metadata": {},
   "source": [
    "And now set some variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "specified-experiment",
   "metadata": {},
   "outputs": [],
   "source": [
    "db_name = '../ispaq_example.db'\n",
    "metric = 'pdf'\n",
    "startDate = '2020-10-01T00:00:00.000000Z'       # Full time is important for retrieving PDFs\n",
    "endDate = '2020-10-15T23:59:59.000000Z'\n",
    "target = 'IU.ANMO.00.BH1.M'\n",
    "\n",
    "startdate = startDate.split('T')[0]\n",
    "enddate = endDate.split('T')[0]\n",
    "filename = f'example3_{target}_{startdate}_{enddate}_PDF.png'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bizarre-yellow",
   "metadata": {},
   "source": [
    "The first step is to create a query that will be used to retrieve the PDFs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "following-defendant",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "SQLcommand = f\"SELECT * FROM {metric} WHERE start = '{startDate}' \" \\\n",
    "             f\"and end = '{endDate}' and (target = '{target}');\"\n",
    "print(SQLcommand)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "internal-durham",
   "metadata": {},
   "source": [
    "Create a connection to the database and run the query, loading it into a pandas dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "twelve-insertion",
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    conn = sqlite3.connect(db_name)\n",
    "    DF = pd.read_sql_query(SQLcommand, conn, parse_dates=['start','end'])\n",
    "    conn.close\n",
    "except:\n",
    "    print(f\"Unable to connect to or find the {metric} table in the database {db_name}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ceramic-exploration",
   "metadata": {},
   "outputs": [],
   "source": [
    "if DF.empty:\n",
    "    print(\"Empty return: there are no PDFs that were retrieved\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "accessory-comparative",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(DF)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "protected-bailey",
   "metadata": {},
   "source": [
    "Sum up the total number of hits for each frequency:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "accompanied-services",
   "metadata": {},
   "outputs": [],
   "source": [
    "for frequency in DF['frequency'].unique():\n",
    "        # Sum hits for total column\n",
    "        DF.loc[DF['frequency'] == frequency, 'total'] = sum(DF[DF['frequency'] == frequency]['hits'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "secret-civilization",
   "metadata": {},
   "source": [
    "For each frequency-power bin, calculate what percentage of the total hits for that frequency are at that power. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "outside-january",
   "metadata": {},
   "outputs": [],
   "source": [
    "DF['percent'] = DF['hits'] / DF['total'] * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tracked-click",
   "metadata": {},
   "source": [
    "Create a minimum range of powers (Y-axis) for better viewing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "optimum-community",
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 = int(min(DF['power'].unique()))\n",
    "p2 = int(max(DF['power'].unique()))\n",
    "if p1 > -190:\n",
    "    p1 = -190\n",
    "if p2 < -90:\n",
    "    p2 = -90\n",
    "    \n",
    "powers = sorted(range(p1,p2+1), reverse=True)\n",
    "freqs = sorted(DF['frequency'].unique(),reverse = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "conceptual-diabetes",
   "metadata": {},
   "source": [
    "Create a new dataframe for plotting: rows are powers, columns are periods, value is percent of hits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "banned-spring",
   "metadata": {},
   "outputs": [],
   "source": [
    "plotDF = pd.DataFrame(0,index=powers,columns=freqs)\n",
    "nonZeroFreqs=[]\n",
    "for power in powers:\n",
    "    for freq in freqs:\n",
    "        value = DF[(DF['frequency']==freq) & (DF['power']== power)]['percent'].values\n",
    "        try:\n",
    "            plotDF.loc[power,freq] = value[0]\n",
    "            if value[0] != 0:\n",
    "                # Keep track of the frequencies that have hits, for axes limits\n",
    "                nonZeroFreqs.append(freq)\n",
    "        except:\n",
    "            continue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "revised-workstation",
   "metadata": {},
   "source": [
    "Matplotlib imshow takes a list (matrix) of values, so convert the dataframe to a list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "grave-antique",
   "metadata": {},
   "outputs": [],
   "source": [
    "plotList = plotDF.values.tolist()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "opposed-traveler",
   "metadata": {},
   "source": [
    "And now we set up some plotting options:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "spectacular-incident",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up plotting -- color map\n",
    "cmap = plt.get_cmap('gist_heat', 3000)  # You can change the colormap here\n",
    "cmaplist = [cmap(i) for i in range(cmap.N)]\n",
    "\n",
    "# convert the first nchange to fade from white, so that anywhere without any hits (or very few) is white\n",
    "nchange = 100\n",
    "for i in range(nchange):\n",
    "    first = cmaplist[nchange][0]\n",
    "    second = cmaplist[nchange][1]\n",
    "    third = cmaplist[nchange][2]\n",
    "    scaleFactor = (nchange-1-i)/float(nchange)\n",
    "\n",
    "    df = ((1-first) * scaleFactor) + first\n",
    "    ds = ((1-second)* scaleFactor) + second\n",
    "    dt = ((1-third) * scaleFactor) + third\n",
    "\n",
    "    cmaplist[i] = (df, ds, dt, 1)\n",
    "\n",
    "cmaplist[0] = (1,1,1,1)\n",
    "cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "photographic-angle",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up plotting -- axis labeling and ticks\n",
    "periodPoints = [0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000]\n",
    "freqPoints = [1/float(i) for i in periodPoints]\n",
    "xfilter = [(i <= freqs[0]) and (i >= freqs[-1]) for i in freqPoints]\n",
    "xlabels = [i for (i, v) in zip(freqPoints, xfilter) if v]\n",
    "xticks = [find_nearest(freqs, i) for i in xlabels]\n",
    "xlabels = [int(1/i)  if i<=1 else 1/i for i in xlabels]     #convert to period, use decimal only if <1s\n",
    "\n",
    "yticks = [powers.index(i) for i in list(filter(lambda x: (x % 10 == 0), powers))]\n",
    "ylabels = [powers[i] for i in yticks]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "considered-sister",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up plotting -- plot\n",
    "height = ylabels[0] - ylabels[-1]\n",
    "plt.figure(figsize=( 12, (.055*height + .5) ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "express-fancy",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(plotList, cmap=cmap,  vmin=0, vmax=30, aspect=.4, interpolation='bilinear')\n",
    "# Adjust grids, labels, limits, titles, etc\n",
    "plt.grid(linestyle=':', linewidth=1)\n",
    "plt.xlabel('Period (s)',size=18)\n",
    "plt.ylabel(r'Power [$10log_{10}(\\frac{m^2/s^4}{hz}$)][dB]',size=18)\n",
    "\n",
    "plt.xticks(xticks[::-1], xlabels[::-1],size=15)\n",
    "plt.yticks(yticks,ylabels,size=15)\n",
    "\n",
    "xmin=freqs.index(min(nonZeroFreqs))\n",
    "xmax=freqs.index(max(nonZeroFreqs))\n",
    "plt.xlim(xmax,xmin)\n",
    "plt.ylim(max(yticks)+5,min(yticks)-5)\n",
    "\n",
    "plt.title(f\"{target}\\n{startdate} through {enddate}\", size=18)\n",
    "\n",
    "# User has option to include colorbar and/or legend\n",
    "cb = plt.colorbar(fraction=.02)\n",
    "cb.set_label('percent probability',labelpad=5)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "occupied-hands",
   "metadata": {},
   "source": [
    "Save the figure for later use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "behind-asbestos",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.tight_layout()\n",
    "plt.savefig(filename)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
